WGCNA.SFT <- function(datExpr, Title, GeneNumCut){
  datExpr <- datExpr
  type = "unsigned"
  corType = "pearson"
  corFnc = ifelse(corType=="pearson", cor, bicor)
  maxPOutliers = ifelse(corType=="pearson",1,0.05)
  robustY = ifelse(corType=="pearson",T,F)
  m.mad <- apply(datExpr,1,mad)
  datExprVar <- datExpr[which(m.mad > 
                                max(quantile(m.mad, probs=seq(0, 1, GeneNumCut))[2],0.01)),]
  dim(datExprVar)
  datExpr <- as.data.frame(t(datExprVar))
  ## 检测缺失值
  gsg = goodSamplesGenes(datExpr, verbose = 3)
  if (!gsg$allOK){
    # Optionally, print the gene and sample names that were removed:
    if (sum(!gsg$goodGenes)>0) 
      printFlush(paste("Removing genes:", 
                       paste(names(dataExpr)[!gsg$goodGenes], collapse = ",")));
    if (sum(!gsg$goodSamples)>0) 
      printFlush(paste("Removing samples:", 
                       paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));
    # Remove the offending genes and samples from the data:
    dataExpr = dataExpr[gsg$goodSamples, gsg$goodGenes]
  }
  ## sample cluster based on expression values
  nGenes = ncol(datExpr)
  nSamples = nrow(datExpr)
  assign("nGenes",value = nGenes, envir = globalenv())
  assign("nSamples",value = nSamples, envir = globalenv())
  dim(datExpr)
  #head(datExpr)[,1:8]
  ## trait-sample tree
  sampleTree = hclust(dist(datExpr), method = "average")

  pdf(file = paste(Title,"Sample_clustering.pdf",sep = "."),width = 28,height = 5)
  plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")
  dev.off();
 # sample_colors <- numbers2colors(as.numeric(factor(dataTraits$Group)),
 #                               colors = c("green","yellow","red"),
 #                               signed = FALSE)
 # pdf(file = paste(Title,"Trait_Sample_clustering.pdf",sep = "."),width = 28,height = 5)
 # par(mar = c(1,4,3,1),cex=0.8)
 # {plotDendroAndColors(sampleTree, sample_colors,
 #                   groupLabels = colnames(sample),
 #                   cex.dendroLabels = 0.8,
 #                   marAll = c(1, 4, 3, 1),
 #                   cex.rowText = 0.01,
 #                   main = "sample dedrogram and trait heatmap")
 # }
 # dev.off();
  # export tree as nwk
  mytree <- as.phylo(sampleTree)
  write.tree(mytree,file = "samplecluster.nwk")
  # sft plot
  powers = c(c(1:10), seq(from = 12, to=30, by=2))
  sft = pickSoftThreshold(datExpr, powerVector=powers, 
                          networkType=type, verbose=5)
  
  pdf(file = paste(Title,"SFTPlot.pdf",sep = "."),width = 10,height = 7)
  par(mfrow = c(1,2))
  cex1 = 0.9
  # 横轴是Soft threshold (power)，纵轴是无标度网络的评估参数，数值越高，
  # 网络越符合无标度特征 (non-scale)
  {plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
         xlab="Soft Threshold (power)",
         ylab="Scale Free Topology Model Fit,signed R^2",type="n",
         main = paste("Scale independence"))
    text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
         labels=powers,cex=cex1,col="red")
    # 筛选标准。R-square=0.85
    abline(h=0.90,col="red")
    abline(h=0.85,col="green")
    # Soft threshold与平均连通性
    plot(sft$fitIndices[,1], sft$fitIndices[,5],
         xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
         main = paste("Mean connectivity"))
    text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, 
         cex=cex1, col="red")
  }
  dev.off();
  power = sft$powerEstimate
  power
  assign("datExpr",value = datExpr, envir = globalenv())
  assign("power",value = power, envir = globalenv())
  assign("sampleTree", value = sampleTree, envir = globalenv())
}
